setwd("C:\\Users\\aniverb\\Documents\\Grad_School\\JHU\\700")

tv=read.table("trades_pivot_2016-11-18.txt", sep="\t", header=T) 
pv=read.table("trades_pivot_vol_2016-11-18.txt", sep="\t", header = T)
trade_data=read.table("trades_count_regression_2016-11-18.txt", sep="\t", header=T)
tv_pv_comb=cbind(tv[,1:16], pv[,4:16], tv[17], pv[17])
tv_pv_comb["Day"]=trade_data["Day"]
head(tv_pv_comb)
##   ProductName   Maturity       Date    X6    X7    X8    X9   X10   X11
## 1          6E 2016-06-13 2016-05-02  4743  4371 14131 19188 26520 10183
## 2          6E 2016-06-13 2016-05-03  7761  9281 23010 18470 27603 18569
## 3          6E 2016-06-13 2016-05-04  4065 10092 31111 16142 41140 12268
## 4          6E 2016-06-13 2016-05-05 11257  9318 17680 22189 30305 12814
## 5          6E 2016-06-13 2016-05-06  3582  4173 58424 29355 16648  7940
## 6          6E 2016-06-13 2016-05-09  4806  9660 12547  9459 14289  6987
##    X12   X13   X14  X15  X16 X17 X18  X6.1  X7.1   X8.1   X9.1  X10.1
## 1 9504  5228  6112 2577 2616   0   0 4.195 2.388  5.414 15.202 10.399
## 2 8871 10705  7300 3884 2149  34   9 2.845 4.866  5.204  6.455 12.777
## 3 4785  4985  6424 2846 1784  15   0 2.709 5.932 11.152  7.811  9.139
## 4 7126  6431  4706 3802 1169   1   0 3.966 3.955  4.532  8.688  7.442
## 5 5493  6323 10189 4226 1304   1   0 3.432 1.998 19.902  8.070  8.105
## 6    0     0     0    0    0  16   0 2.251 4.298  6.758  4.963  7.391
##   X11.1 X12.1 X13.1 X14.1 X15.1 X16.1 X17.1 X18.1 DayTradeTotal
## 1 2.440 6.733 3.743 4.343 1.836 5.190    NA    NA        105173
## 2 3.812 6.492 2.912 6.776 1.797 3.144 0.866 0.837        137646
## 3 4.576 2.625 3.157 3.281 2.045 2.589 0.354    NA        135657
## 4 3.817 4.305 3.403 2.657 1.482 1.047    NA    NA        126798
## 5 2.980 2.001 5.503 4.956 1.861 1.606    NA    NA        147658
## 6 3.536    NA    NA    NA    NA    NA 0.635    NA         57764
##   DayVolatility Day
## 1        16.844   1
## 2        23.703   2
## 3        12.542   3
## 4        14.243   4
## 5        17.925   5
## 6         8.338   6
pvals=c()
prod=c()
products=as.vector(unique(tv_pv_comb$ProductName))
significant=list(c("ZQ", "2016-08-31"), c("YM", "2017-03-17"), c("HO", "2016-11-30"), c("ZF", "2016-12-30"), c("RB", "2016-07-29"))
dfInt=data.frame(Products=c("HO", "6E", "CL", "ES", "KE", "ZL", "CL", "GC", "RB", "GE", "GE", "GE", "UB", "ZB", "ZM", "ZQ"),  Maturities=c("2016-09-30","2017-06-19", "2016-07-20", "2016-09-16", "2016-12-14","2016-12-14", "2017-01-20", "2016-10-27", "2016-07-29", "2016-09-19", "2016-07-18", "2016-11-14", "2016-09-21", "2016-09-21", "2016-10-14", "2016-07-29"))

inDF=function(x, product, maturity){
  if ((x[1]==product) & (x[2]==maturity)){
    return(TRUE)
  }else{
    return(FALSE)
  }
}
for (product in products){ #in results, put 3 plots on 1 line
  tv_pv_combS=subset(tv_pv_comb, ProductName==product)
  maturities=as.vector(unique(tv_pv_combS$Maturity))
  for (maturity in maturities){
    tv_pv_combM=subset(tv_pv_combS, Maturity==maturity)
      if (nrow(tv_pv_combM)>=30){
        cors=apply(tv_pv_combM, 1, function(x){cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use="na.or.complete")})
        dfCors=data.frame(tv_pv_combM[c("Date", "Day", "DayTradeTotal")], cors=as.vector(cors))
        tradeQuant=quantile(dfCors$DayTradeTotal, c(.6, .9)) 
        tv_pv_combMQ=dfCors[(dfCors$DayTradeTotal>tradeQuant[1]) & (dfCors$DayTradeTotal<tradeQuant[2]),] #selecting high volume trades
        mod=lm(tv_pv_combMQ$cors~tv_pv_combMQ$Day) 
        cor_hat=predict(mod, tv_pv_combMQ["Day"])
        tv_pv_combMQ["cor_hat"]=cor_hat
        
        if ((product==significant[[1]][1]) & (maturity==significant[[1]][2]) | (product==significant[[2]][1]) & (maturity==significant[[2]][2]) | (product==significant[[3]][1]) & (maturity==significant[[3]][2])| (product==significant[[4]][1]) & (maturity==significant[[4]][2])| (product==significant[[5]][1]) & (maturity==significant[[5]][2])){
          png(paste("Coeff_fit", product, maturity, ".png", sep=""), width = 6, height = 4, units = 'in', res = 400)
          plot(tv_pv_combMQ$Day, tv_pv_combMQ$cors, ylim=c(-1,1), xlab="Day", ylab="correlation coefficient", main=paste(product, "Maturity:", maturity), xaxt="n")
          lines(tv_pv_combMQ$Day, tv_pv_combMQ$cor_hat, col="red")
          at=seq(min(tv_pv_combMQ$Day), max(tv_pv_combMQ$Day), 1)
          labs=rep(NA, length(at))
          labs[which(at%in%tv_pv_combMQ$Day)]=as.vector(tv_pv_combMQ$Date)
          axis(1, at=at, labels=labs, tick=F)
          dev.off()
        }
        
        plot(tv_pv_combMQ$Day, tv_pv_combMQ$cors, ylim=c(-1,1), xlab="Day", ylab="correlation coefficient", main=paste(product, "Maturity:", maturity), xaxt="n")  
        lines(tv_pv_combMQ$Day, tv_pv_combMQ$cor_hat, col="red")
        at=seq(min(tv_pv_combMQ$Day), max(tv_pv_combMQ$Day), 1)
        labs=rep(NA, length(at))
        labs[which(at%in%tv_pv_combMQ$Day)]=as.vector(tv_pv_combMQ$Date)
        axis(1, at=at, labels=labs, tick=F)
        
        stat=summary(mod)
        stat=summary(mod)
        cof=stat$coefficients
        pval=cof[,4][2]
        pvals=append(pvals, pval, after=length(pvals))
        prod=append(prod, paste(product, maturity), after=length(prod))
        
        trades=as.vector(as.matrix(t(tv_pv_combM[,4:16])))
        priceVol=as.vector(as.matrix(t(tv_pv_combM[,17:29])))
        df=data.frame(trades, priceVol)[order(trades),]
        df=na.omit(df)
        n=nrow(df)
        dfTrim=df[ceiling(.05*n):round(.95*n),]
        dfTrimO=dfTrim[order(dfTrim$priceVol),]
        n2=nrow(dfTrimO)
        dfTrimP=dfTrimO[ceiling(.05*n2):round(.95*n2),]
        
        plot(dfTrimP, xlab="trade volume hourly", ylab="price volatility hourly", main=paste(product, "Maturity:", maturity), pch=20) 
        
        daily=tv_pv_combM[,30:31][order(tv_pv_combM$DayTradeTotal),]
        daily=na.omit(daily)
        nD=nrow(daily)
        dailyTrim=daily[ceiling(.05*nD):round(.95*nD),]
        dailyTrimO=dailyTrim[order(dailyTrim$DayVolatility),]
        nD2=nrow(dailyTrimO)
        dailyTrimP=dailyTrimO[ceiling(.05*nD2):round(.95*nD2),]
        plot(dailyTrimP, xlab="trade volume daily", ylab="price volatility daily", main=paste(product, "Maturity:", maturity), pch=20)
        
        t=apply(dfInt, 1, inDF, product, maturity)
        if (TRUE %in% t){
          png(paste(product, maturity, "hourly.png", sep=""), width = 6, height = 4, units = 'in', res = 300)
          plot(dfTrimP, xlab="trade volume hourly", ylab="price volatility hourly", main=paste(product, "Maturity:", maturity), col="pink3", pch=20)
          dev.off()
          png(paste(product, maturity, "daily.png", sep=""), width = 6, height = 4, units = 'in', res = 300)
          plot(dailyTrimP, xlab="trade volume daily", ylab="price volatility daily", main=paste(product, "Maturity:", maturity), col="pink3", pch=20)
          dev.off()
        }
    }
  }
}

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

## Warning in cor(as.numeric(x[4:16]), as.numeric(x[17:29]), use =
## "na.or.complete"): the standard deviation is zero

pvalsO=order(pvals)
corPvals=data.frame(prod, pvals)[pvalsO,]
nPvals=nrow(corPvals)
nPvals
## [1] 126
corPvals["id"]=1:nPvals #B-H Procedure
corPvals["q"]=nPvals*corPvals$pvals/corPvals$id
f=c()
for(i in 1:nPvals){f=append(f, min(corPvals$q[i:nPvals]), length(f))}
corPvals["f"]=f
corSig=corPvals[which(corPvals$pvals<.05),]
nSig=nrow(corSig)
nSig
## [1] 16
corSig #5 rejected at FDR=.20 
##              prod       pvals id         q         f
## 111 ZQ 2016-08-31 0.001326200  1 0.1671012 0.1181514
## 78  YM 2017-03-17 0.002226816  2 0.1402894 0.1181514
## 36  HO 2016-11-30 0.002813129  3 0.1181514 0.1181514
## 90  ZF 2016-12-30 0.004494289  4 0.1415701 0.1415701
## 59  RB 2016-07-29 0.006932494  5 0.1746988 0.1746988
## 43  KE 2017-05-12 0.014340904  6 0.3011590 0.3011590
## 6   CL 2016-07-20 0.018004212  7 0.3240758 0.3240758
## 2   6E 2016-12-19 0.021064406  8 0.3317644 0.3302329
## 103 ZM 2016-12-14 0.023588068  9 0.3302329 0.3302329
## 32  HO 2016-07-29 0.027519890 10 0.3467506 0.3467506
## 24  GC 2017-04-26 0.032399027 11 0.3711161 0.3546173
## 29  GE 2016-12-19 0.034721438 12 0.3645751 0.3546173
## 53  NG 2017-04-26 0.039869811 13 0.3864305 0.3546173
## 3   6E 2017-03-13 0.040719034 14 0.3664713 0.3546173
## 115 ZQ 2016-12-30 0.042216346 15 0.3546173 0.3546173
## 55  NQ 2016-09-16 0.047810841 16 0.3765104 0.3765104
#trimming seems to not have worked for ZQ and RB since they are close to maturity

set.seed(3-3-17)
u=sort(runif(length(pvals)))
#png("unifQQ.png")
plot(u, corPvals$pvals, ylab="p-values", xlab="uniform", main="Q-Q Plot")
abline(0,1, col="red")

#dev.off()